Warning: package 'rstan' was built under R version 4.1.2
Loading required package: StanHeaders
Loading required package: ggplot2
Warning: package 'ggplot2' was built under R version 4.1.2
rstan (Version 2.21.5, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Show code
library(posterior)
Warning: package 'posterior' was built under R version 4.1.2
This is posterior version 1.2.1
Attaching package: 'posterior'
The following objects are masked from 'package:rstan':
ess_bulk, ess_tail
The following objects are masked from 'package:stats':
mad, sd, var
Show code
library(ggplot2)library(mgcv)
Loading required package: nlme
This is mgcv 1.8-35. For overview type 'help("mgcv-package")'.
Show code
library(dimreduce)
This is dimreduce version 0.2.2
Show code
source("R/rstan.R")source("R/search.R")source("R/project.R")source("R/project_old.R") # for project_sigma()source("R/simulate.R")
Warning: package 'MASS' was built under R version 4.1.2
# Run sim and create data setupdat <-simulate(n, rho, rel_true, sigma)n <-length(dat$y)splt <-create_split(n, 0.5)ds <-list(dat = dat, split = splt)
Fitting reference model
Show code
# Model setup# Create the Stan modelmod <- rstan::stan_model("stan/model_sample.stan")ms <-list(stan_model = mod,B =24,scale_bf =1.5)dat_df <-create_data_frame(ds, test =FALSE)fmf1 <-fit_full_model(ds, ms,chains = CHAINS, iter = ITER, control = CN,thin =10,seed =507949442)
Warning: The largest R-hat is 1.45, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat
Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess
Sampling done.
Show code
# Use SPCA reference model# fmf2 <- fit_full_model(ds, ms,# chains = CHAINS, iter = ITER, control = CN,# thin = 10, nctot = 5# )
Ranking the variables
Show code
# Rank and print the orderingcv1 <-comp_vars(fmf1$fit, fmf1$J)path <-sort(cv1, index.return = T, decreasing = T)$ixprint(path)
[1] 1 5 2 8 3 7 6 4
Visualizing the reference model posterior
Show code
pa <-plot(fmf1$fit, pars =c("alpha")) +ggtitle("Magnitudes")
ci_level: 0.8 (80% intervals)
outer_level: 0.95 (95% intervals)
Show code
pb <-plot(fmf1$fit, pars =c("ell")) +ggtitle("Lengthscales")
# Run selection with projection predictive methodoptions(pp.threshold =Inf) # no stopping conditionoptions(min.sp =1e-1) # lower bound for smoothing parameter# Use HS GP basisoptions(bs.type ="hs")sel_pp_hs4 <-selection_pp(fmf1, path)
Performing forward search.
Running a ProjectionForwardSearch with option = gam
Stopping threshold for explained variance = Inf
Projecting to submodel []
Step 1/8.
* using pre-defined search path
Projecting to submodel [1]
* 1 - kl_current/kl_empty = 0.5188
Step 2/8.
* using pre-defined search path
Projecting to submodel [1 5]
* 1 - kl_current/kl_empty = 0.8724
Step 3/8.
* using pre-defined search path
Projecting to submodel [1 5 2]
* 1 - kl_current/kl_empty = 0.9732
Step 4/8.
* using pre-defined search path
Projecting to submodel [1 5 2 8]
* 1 - kl_current/kl_empty = 0.9896
Step 5/8.
* using pre-defined search path
Projecting to submodel [1 5 2 8 3]
* 1 - kl_current/kl_empty = 0.9942
Step 6/8.
* using pre-defined search path
Projecting to submodel [1 5 2 8 3 7]
* 1 - kl_current/kl_empty = 0.9966
Step 7/8.
* using pre-defined search path
Projecting to submodel [1 5 2 8 3 7 6]
* 1 - kl_current/kl_empty = 0.9988
Step 8/8.
* using pre-defined search path
Projecting to submodel [1 5 2 8 3 7 6 4]
* 1 - kl_current/kl_empty = 1
Search done.
Show code
# Use thin plate spline basisoptions(bs.type ="tp")sel_pp_tp <-selection_pp(fmf1, path)
Performing forward search.
Running a ProjectionForwardSearch with option = gam
Stopping threshold for explained variance = Inf
Projecting to submodel []
ell_extreme <-extract(fmf1$fit, "ell")$ell[d_idx,]alpha_extreme <-extract(fmf1$fit, "alpha")$alpha[d_idx,]cat("Printing the reference model lengthscales in the most extreme draw\n")
Printing the reference model lengthscales in the most extreme draw
Show code
print(round(ell_extreme, 2))
[1] 1.31 0.87 1.97 0.59 1.19 1.26 2.21 2.16
Show code
cat("Printing the reference model magnitude params in the most extreme draw\n")
Printing the reference model magnitude params in the most extreme draw